#ifndef AnalysisUtils_h
#define AnalysisUtils_h

#include "BaseDataClasses.hxx"
#include "DetectorDefinition.hxx"
#include "SystematicSource.hxx"
#include "BasicUtils.hxx"

/// This namespace contains useful functions for analyses. Functions include
/// those related to fiducial volumes, accessing tracks with specific
/// characteristics, and more.
namespace anaUtils{

    //----- Utility functions -------------
    void ClearArray(AnaTrackB** objArr, int nObj);
    void FillArray(AnaTrackB** tgtArr, AnaTrackB** srcArr, int nObj);
    void CreateArray(AnaTrackB** &tgtArr, int nObj);

    void ClearArray(AnaTrueTrackB** objArr, int nObj);
    void FillArray(AnaTrueTrackB** tgtArr, AnaTrueTrackB** srcArr, int nObj);
    void CreateArray(AnaTrueTrackB** &tgtArr, int nObj);

    /// Get The sofware version from the path of the package
    std::string GetSoftwareVersionFromPath(const std::string& path);

    /// Calculate the distance between two points
    float GetSeparationSquared(const Float_t* pos1, const Float_t* pos2);

    /// Whether the specified position is in the volume of the given detector.
    /// Accepted detectors are kFGD1, kFGD2, kP0D, kDSECal, k(Top, Bottom, Left, Right)T(P)ECal (or SMRD)
    bool InDetVolume(SubDetId::SubDetEnum det, const Float_t* pos);

    /// Whether the specified position is in the fiducial volume of the given detector.
    /// Accepted detectors are kFGD1, kFGD2, kP0D, kDSECal, k(Top, Bottom, Left, Right)T(P)ECal (or SMRD)
    /// The fiducial volume is specified using the FVdefmin and FVdefmax vectors. These
    /// vectors are the amount of space to *exclude*, compared to the nominal side of
    /// the detector.
    bool InFiducialVolume(SubDetId::SubDetEnum det, const Float_t* pos, const Float_t* FVdefmin, const Float_t* FVdefmax);
    bool InFiducialVolume(SubDetId::SubDetEnum det, const Float_t* pos);

    /// Returns the material of FGD2 in which the Z position is.
    // to do: put geometry in DetectorDefinition and implement x,y
    int GetMaterialInFgd2(const Float_t* pos);

    /// Returns the detector in which the position is.
    SubDetId::SubDetEnum GetDetector(const Float_t* pos);

    /// Get all delayed time bins as Michel Electron candidates
    int GetFGDMichelElectrons(const AnaEventB& event, const SubDetId::SubDetEnum det, AnaFgdTimeBinB** arr, bool prod5Cut = 0);

    /// Get the number of true michel electrons
    Int_t GetNMichelElectrons(const AnaTrueVertexB& trueVertex, SubDetId::SubDetEnum det = SubDetId::kFGD1);

    /// Return the true linear length traversed in the TPC
    Float_t GetTrueLinearLengthInTPC(const AnaTrueTrackB& trueTrack, Float_t& distz);

    /// For tracks that start in the FGD, get the closest TPC in the direction of the track.
    SubDetId::SubDetEnum GetClosestTPC(const AnaTrackB& track);

    /// Get the vector of AnaSubTrack segment that uses the specified detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.  
    /// Returns the number of entries in the input array, the number of tracks found
    int GetSegmentsInDet(const AnaTrackB& track, SubDetId::SubDetEnum det, AnaSubTrackB* selTracks[]);

    /// Method to get the subtrack with most nodes in a given detector
    AnaSubTrackB* GetSegmentWithMostNodesInDet(const AnaTrackB& track, SubDetId::SubDetEnum det);

    /// Combined function to address NuMu selection needs as efficiently as possible - gets the TPC segment with the most nodes in the TPC closest to the start of the global track
    AnaSubTrackB* GetSegmentWithMostNodesInClosestTpc(const AnaTrackB& track);

    /// Get the AnaSubTrack segment that uses the specified detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    AnaSubTrackB* GetSegmentInDet(const AnaTrackB& track, SubDetId::SubDetEnum det);

    /// Get the number of AnaSubTrack segments that use the specified detector system.
    Int_t GetNSegmentsInDet(const AnaTrackB& track, SubDetId::SubDetEnum det);

    /// Compare the NHits of two tracks. Return whether t1 is higher
    /// NHts than t2.
    bool CompareNHits(AnaTrackB* t1, AnaTrackB* t2);

    bool HasTrackUsingTPC(const AnaEventB& event);
    bool HasTrackUsingDet(const AnaEventB& event, SubDetId::SubDetEnum det);

    /// Get all the true traj. in the bunch, 
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTrajInBunch(const AnaEventB& event, AnaTrueTrackB* traj[]);
    /// Get all the true traj. in the bunch that are charged 
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllChargedTrajInBunch(const AnaEventB& event, AnaTrueTrackB* traj[]);
    /// Get all the true traj. in the bunch that are charged and crosses the TPC 
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllChargedTrajInTPCInBunch(const AnaEventB& event, AnaTrueTrackB* traj[]);
    /// Get all the true traj. in the bunch that are charged and crosses the FGD
    /// Returns the number of entries in the input array, the number of tracks found
  int GetAllChargedTrajInFGDInBunch(const AnaEventB& event, AnaTrueTrackB* traj[],SubDetId::SubDetEnum det);
    /// Get all the true traj. in the bunch that are charged and crosses the TPC with a length > 1/4 of the TPC 
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllBigEnoughChargedTrajInTPCInBunch(const AnaEventB& event, AnaTrueTrackB* chargedtrajInBunch[]);
    /// Get all the true traj. in the bunch that are charged and crosses the TPC and the FGD (TPC1-FGD1, FGD1-TPC2, TPC2-FGD2, FGD2-TPC3)
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllChargedTrajInTPCFGDInBunch(const AnaEventB& event, AnaTrueTrackB* chargedtrajInBunch[]);
    /// Get all the true traj. in the bunch that are charged and crosses the the FGD but not the TPC
    /// Returns the number of entries in the input array, the number of tracks found
  int GetAllChargedTrajInFGDAndNoTPCInBunch(const AnaEventB& event, AnaTrueTrackB* chargedtrajInBunch[],SubDetId::SubDetEnum det);

    /// Get all the tracks using a specific detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingDet(const AnaEventB& event, SubDetId::SubDetEnum det, AnaTrackB* selTracks[]);

    /// Get all the tracks not using a specific detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksNotUsingDet(const AnaEventB& event, SubDetId::SubDetEnum det, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the P0D, sorted by decreasing momentum.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingP0D(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the FGD, sorted by decreasing momentum.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingFGD(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the TPC, sorted by decreasing momentum.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingTPC(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the TPC or the FGD, sorted by decreasing number of hits
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingFGDorTPC(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the ECAL, sorted by decreasing momentum.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingECAL(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the SMRD, sorted by decreasing momentum.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingSMRD(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Access function to get all the tracks in the bunch that use the FGD and no TPC, sorted by decreasing NHits.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingFGDAndNoTPC(const AnaEventB& event, AnaTrackB* selTracks[]);

    /// Returns the number of tracks using both the TPC and the subdetector 'det'.
    int GetNTracksUsingTPCAndDet(const AnaEventB& event, SubDetId::SubDetEnum det);

    /// Whether a track uses all the specified detector systems.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    bool TrackUsesDets(const AnaTrackB& track, SubDetId::SubDetEnum dets[], int nDets);

    /// Whether a track uses the specified detector system.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    bool TrackUsesDet(const AnaTrackB& track, SubDetId::SubDetEnum det);

    /// Whether a subtrack uses the specified detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    bool TrackUsesDet(const AnaSubTrackB& track, SubDetId::SubDetEnum det);

    /// Get all the tracks using only a specific detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetAllTracksUsingOnlyDet(const AnaEventB& event, SubDetId::SubDetEnum det, AnaTrackB* selTracks[]);

    /// Whether a track uses only the specified detector system.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    bool TrackUsesOnlyDet(const AnaTrackB& track, SubDetId::SubDetEnum det);

    /// Whether a subtrack uses only the specified detector.
    /// See SubDetId::SubDetEnum for the detector numbering convention.
    bool TrackUsesOnlyDet(const AnaSubTrackB& track, SubDetId::SubDetEnum det);

    /// Merge the two vectors of AnaTrackBs into a single vector, with no duplicated tracks.
    std::vector<AnaTrackB*> MergeUniqueTracks(const std::vector<AnaTrackB*>& tracks1, const std::vector<AnaTrackB*>& tracks2);

    /// Return a scale factor for the magnet intensity, based on the ND280 run number.
    /// THIS IS CURRENTLY DISABLED AND WILL ALWAYS RETURN 1!
    Float_t MagnetIntensityCorrection(Int_t Run);

    /// Function to recompute the pull for a TPC track segment
    Float_t ComputeTPCPull(const AnaTpcTrackB &track, const std::string& particle);

    /// Function to recompute the pull for a TPC track segment for all hypotheses
    void ComputeTPCPull(const AnaTpcTrackB &track, Float_t* pulls);

    /// Function to recompute all the pull for a TPC track segment and save them into the segment
    void RecomputeTPCPulls(AnaTpcTrackB &track);

    /// Compute the expected TPC dEdx for a track, using a different
    /// parameterisation than the one in tpcRecon.
    Float_t ExpectedTPCdEdx(const AnaTpcTrackB &track, const std::string& particle);

    /// The likelihood of a track being a given particle hypothesis, based on the
    /// pull values of all TPC segments in the track.
    ///
    /// hypo is one of:
    /// * 0 - muon
    /// * 1 - electron
    /// * 2 - proton
    /// * 3 - pion
    Float_t GetPIDLikelihood(const AnaTrackB& track, Int_t hypo, bool prod5Cut = 0);

    /// Get all likelihoods
    void GetPIDLikelihood(const AnaTrackB&, Float_t* hypo, bool prod5Cut = 0);

    /// Get the likelihood for MIP: (like_mu+like_pi)/(1-like_p)
    Float_t GetPIDLikelihoodMIP(const AnaTrackB &track);

    /// A function that is not currently used, and will be documented when it is.
    Float_t GetPIDPrior(const AnaTrackB& track, Int_t hypo);

    /// compute the inverse transverse momentum given an AnaTrackB
    Float_t ComputeInversePT(const AnaTrackB &track);

    /// compute the inverse transverse momentum given an AnaTrackB and the momentum
    Float_t ComputeInversePT(const AnaTrackB &track, Float_t mom);

    /// compute the inverse transverse momentum given an AnaTpcTrackB
    Float_t ComputeInversePT(const AnaTpcTrackB &track);

    /// compute the inverse transverse momentum for a true track
    Float_t ComputeInversePT(const AnaTrueTrackB &track);

    /// compute the total momentum given the track and the inverse transverse momentum
    Float_t ComputeMomentumFromInversePT(const AnaTrackB &track, Float_t PTinv);

    /// compute the total momentum given the track and the inverse transverse momentum
    Float_t ComputeMomentumFromInversePT(const AnaTpcTrackB &track, Float_t PTinv);

    /// Compute the reconstructed neutrino energy assuming CCQE kinematics. Uses the "basic"
    /// formula which doesn't consider neutron/proton mass difference or binding energy.
    ///
    /// Formula is E_rec = 0.5 * (2 * M_p * E_l - m_l^2) / (M_p - E_l + p_l * costheta_l),
    /// where M_ is mass, E_ is energy, p_ is momentum, p is proton, l is lepton.
    Float_t NeutrinoERecCCQE(Float_t mom_lepton, Float_t mass_lepton, Float_t costheta_lepton);


    /// Compute the reconstructed neutrino energy assuming CCQE kinematics. 
    /// It considers neutron/proton mass difference and binding energy.
    Float_t ComputeRecNuEnergyCCQE(Float_t mom_lepton, Float_t mass_lepton, Float_t costheta_lepton);

    /// Compute the reconstructed Q2 assuming CCQE kinematics to compute neutrino energy. 
    /// It considers neutron/proton mass difference and binding energy.
    Float_t ComputeRecQ2CCQE(Float_t mom_lepton, Float_t mass_lepton, Float_t costheta_lepton);

    /// Compute the Q2 using muon variables and neutrino energy
    /// It considers neutron/proton mass difference and binding energy.
    Float_t ComputeQ2(Float_t mom_lepton, Float_t mass_lepton, Float_t costheta_lepton, Float_t energy_neutrino);

    /// Get the AnaTrueTrackB in the current spill with the given ID. Return NULL if it can't be found.
    AnaTrueTrackB* GetTrueTrackByID(const AnaEventB& event, int ID);

    /// Whether this track has an ecal segment in the TrackerECal or DsECal.
    bool UsesTrackerDsEcal(AnaTrackB* track);

    /// Whether this ecal segment is in the TrackerECal or DsECal (i.e. not in the P0DECal).
    bool UsesTrackerDsEcal(AnaEcalTrackB* ecal);

    /// Get the Tracker ECal and DsECal segments from the given track.
    /// Returns the number of entries in the input array, the number of tracks found
    int GetTrackerDsEcals(AnaTrackB* track, AnaEcalTrackB* selTracks[]);

    /// Get the most dowstream TPC segment of the track.
    AnaTpcTrackB* GetTPCBackSegment(const AnaTrackB* track);

    /// Get the most upstream ECal component of a track, based on either the start or
    /// end position. Fill the pos variable with the position used. Return NULL if not found.
    AnaEcalTrackB* GetMostUpstreamECalSegment(AnaTrackB* track);

    /// Whether the ECal object is "contained" within the ECal volume. This is
    /// intended for deciding whether to trust the ECal EM energy or not.
    /// At the moment this is based on the start position of the shower. It will
    /// be updated to a more rigorous definition when production 5F files have been
    /// studied more closely.
    //  bool IsEcalContained(AnaEcalTrackEcalPID* EcalSegment);

    /// Returns the run period (sequentially: 0,1,2,3,4,5 ...)
    int GetRunPeriod(int run);
  
  /// Returns the sans muon modes (sequentially: 0,1,2)
    int GetSandMode(int run);
  
  /// Return true if we are analysing fgd1/fgd2...
  bool IsFGD1Analysis(const AnaEventB& event);
  bool IsFGD2Analysis(const AnaEventB& event);

    /// Return an integer corresponding to the array index of the track in the old local detector enumeration
    int GetLocalDetEnum(SubDetId::SubDetEnum det, int i);

    /// Compute all event weights (FluxWeight, etc. not the systematic ones !!!) and return the vector
    std::vector<Float_t> ComputeWeights(const AnaEventB& event);

    /// Initialise the event systematics box
  void InitializeEventSystematics(const AnaEventB& event, bool initialisePionSystematics=true,bool doFGD1=true,bool doFGD2=false);

    /// Method to load external ipon cross section data into arrays
    template<int nIntBins, int nMomBins>
    void LoadPionCrossSections(char* inputFileName, Float_t (&xsec_array)[nIntBins][nMomBins], Float_t (&err_array)[nIntBins][nMomBins]){
        ifstream inputFile(inputFileName, std::ios::in);

        if (!inputFile) {
            std::cerr << "Cannot open input file '" << inputFileName << "'. Exit!" << std::endl;
            exit(0);
        }

        std::string ctype, cintBin, cvalue;
        char* pEnd;
        while (inputFile >> ctype >> cintBin) {
            int intBin = strtod(cintBin.c_str(),&pEnd);
            int type = strtod(ctype.c_str(),&pEnd);

            if(intBin >= nIntBins) std::cout << "Error - interaction bin larger than input array size" << std::endl;

            for (int momBin = 0;momBin<nMomBins;momBin++){
                inputFile >> cvalue;
                Float_t value = strtod(cvalue.c_str(),&pEnd);
                if (type==0){
                    xsec_array[intBin][momBin] =  value;
                }
                else
                    err_array[intBin][momBin] =  value;
            }
        }

        inputFile.close();
    }

}
#endif

//  LocalWords:  ifndef
